* ********************************************************************************************
* The Effects of Weather Shocks on Economic Activity: What are the Channels of Impact?
* Sebastian Acevedo, Mico Mrkaic, Natalija Novta, Evgenia Pugacheva, Petia Topalova
* With support from Gavin Asdorian, Olivia Ma, Jilun Xing and Yuan Zeng 
* Replication files for Table 1, column 8
* ********************************************************************************************
clear all
set more off
set matsize 11000

cd "C:\Users\..\replication" // set your working directory to where the replication folder is saved

* the Y variable that will be used in regressions: ngdprpc
global Y ngdprpc

* Horizons
global k 7

* ****************************************************************************************
* Our data - create extra variables
* ****************************************************************************************
use "data/country_dataset.dta", clear
drop if year > 2016

gen sqyear = year ^ 2

* create variable for region-year (ry), encode it to be numeric, to be used as fixed effects
gen tempvar1 = wdi_region + " " + string(year) if wdi_region != ""
encode tempvar1, generate(ry)
drop tempvar1

* create country-specific time trend (_IifsXyea*) and time trend squared (_IifsXsqy*)
* dummies for region-year(_Iry*), country (_Iifscode*) and year (_Iyear*)
xi i.ifscode*year i.ifscode*sqyear i.ry i.year

* calculate GDP log difference
gen ln_${Y} = ln(${Y})
gen ${Y}_g = ln_${Y} - l.ln_${Y}

* generate squared terms of temperature and of precipitation
foreach X in CRU_pw50 {
	gen tmp_`X'2 = tmp_`X' ^ 2
	gen precip_`X'2 = precip_`X' ^ 2
}

* ****************************************************************************************
* Variables where to save results
* ****************************************************************************************

gen period = _n-2

foreach X in t t2 p p2 t_ae t_em t_lic p_ae p_em p_lic {
	gen `X'_coef = 0 in 1
	gen `X'_se = 0 in 1
	gen `X'_pvalue = 0 in 1
}

gen nobs = 0 in 1
gen ncountry = 0 in 1
gen r2 = 0 in 1

* ****************************************************************************************
* Regression
* ****************************************************************************************
reg ${Y}_g l(0/$k).tmp_CRU_pw50 l(0/$k).tmp_CRU_pw502  l(0/$k).precip_CRU_pw50 l(0/$k).precip_CRU_pw502 _Iifscode* _Iry*, cluster(ifscode)

* ****************************************************************************************
* Store output
* ****************************************************************************************
replace t_coef = _b[tmp_CRU_pw50] in 2
replace t_se = _se[tmp_CRU_pw50] in 2
replace t_pvalue = 2 * ttail(e(df_r),abs(_b[tmp_CRU_pw50]/_se[tmp_CRU_pw50])) in 2

replace t2_coef = _b[tmp_CRU_pw502] in 2
replace t2_se = _se[tmp_CRU_pw502] in 2
replace t2_pvalue = 2 * ttail(e(df_r),abs(_b[tmp_CRU_pw502]/_se[tmp_CRU_pw502])) in 2

replace p_coef = _b[precip_CRU_pw50] in 2
replace p_se = _se[precip_CRU_pw50] in 2
replace p_pvalue = 2 * ttail(e(df_r),abs(_b[precip_CRU_pw50]/_se[precip_CRU_pw50])) in 2

replace p2_coef = _b[precip_CRU_pw502] in 2
replace p2_se = _se[precip_CRU_pw502] in 2
replace p2_pvalue = 2 * ttail(e(df_r),abs(_b[precip_CRU_pw502]/_se[precip_CRU_pw502])) in 2

forval t = 1/$k {
	replace t_coef = _b[l`t'.tmp_CRU_pw50] in `=`t'+2'
	replace t_se = _se[l`t'.tmp_CRU_pw50] in `=`t'+2'
	replace t_pvalue = 2 * ttail(e(df_r),abs(_b[l`t'.tmp_CRU_pw50]/_se[l`t'.tmp_CRU_pw50])) in `=`t'+2'

	replace t2_coef = _b[l`t'.tmp_CRU_pw502] in `=`t'+2'
	replace t2_se = _se[l`t'.tmp_CRU_pw502] in `=`t'+2'
	replace t2_pvalue = 2 * ttail(e(df_r),abs(_b[l`t'.tmp_CRU_pw502]/_se[l`t'.tmp_CRU_pw502])) in `=`t'+2'	
	
	replace p_coef = _b[l`t'.precip_CRU_pw50] in `=`t'+2'
	replace p_se = _se[l`t'.precip_CRU_pw50] in `=`t'+2'
	replace p_pvalue = 2 * ttail(e(df_r),abs(_b[l`t'.precip_CRU_pw50]/_se[l`t'.precip_CRU_pw50])) in `=`t'+2'

	replace p2_coef = _b[l`t'.precip_CRU_pw502] in `=`t'+2'
	replace p2_se = _se[l`t'.precip_CRU_pw502] in `=`t'+2'
	replace p2_pvalue = 2 * ttail(e(df_r),abs(_b[l`t'.precip_CRU_pw502]/_se[l`t'.precip_CRU_pw502])) in `=`t'+2'	
}


replace nobs = `e(N)'
replace ncountry = `e(N_clust)'
replace r2 = `e(r2_a)'

*** Temperature
* Lincom AE
local for_lincom _b[tmp_CRU_pw50] + (2 * _b[tmp_CRU_pw502] * 11)
forval t = 0/$k {
	lincom `for_lincom'
	replace t_ae_coef = r(estimate) in `=`t'+2'
	replace t_ae_se = r(se) in `=`t'+2'
	replace t_ae_pvalue = 2 * ttail(r(df),abs(r(estimate)/r(se))) in `=`t'+2'
	
	local for_lincom `for_lincom' + _b[l`=`t'+1'.tmp_CRU_pw50] + (2 * _b[l`=`t'+1'.tmp_CRU_pw502] * 11)
}

* Lincom EM
local for_lincom _b[tmp_CRU_pw50] + (2 * _b[tmp_CRU_pw502] * 22)
forval t = 0/$k {
	lincom `for_lincom'
	replace t_em_coef = r(estimate) in `=`t'+2'
	replace t_em_se = r(se) in `=`t'+2'
	replace t_em_pvalue = 2 * ttail(r(df),abs(r(estimate)/r(se))) in `=`t'+2'
	
	local for_lincom `for_lincom' + _b[l`=`t'+1'.tmp_CRU_pw50] + (2 * _b[l`=`t'+1'.tmp_CRU_pw502] * 22)
}

* Lincom LIC
local for_lincom _b[tmp_CRU_pw50] + (2 * _b[tmp_CRU_pw502] * 25)
forval t = 0/$k {
	lincom `for_lincom'
	replace t_lic_coef = r(estimate) in `=`t'+2'
	replace t_lic_se = r(se) in `=`t'+2'
	replace t_lic_pvalue = 2 * ttail(r(df),abs(r(estimate)/r(se))) in `=`t'+2'
	
	local for_lincom `for_lincom' + _b[l`=`t'+1'.tmp_CRU_pw50] + (2 * _b[l`=`t'+1'.tmp_CRU_pw502] * 25)
}



*** Precipitation
* Lincom AE
local for_lincom _b[precip_CRU_pw50] + (2 * _b[precip_CRU_pw502] * 8)
forval t = 0/$k {
	lincom `for_lincom'
	replace p_ae_coef = r(estimate) in `=`t'+2'
	replace p_ae_se = r(se) in `=`t'+2'
	replace p_ae_pvalue = 2 * ttail(r(df),abs(r(estimate)/r(se))) in `=`t'+2'
	
	local for_lincom `for_lincom' + _b[l`=`t'+1'.precip_CRU_pw50] + (2 * _b[l`=`t'+1'.precip_CRU_pw502] * 8)
}

* Lincom EM
local for_lincom _b[precip_CRU_pw50] + (2 * _b[precip_CRU_pw502] * 9)
forval t = 0/$k {
	lincom `for_lincom'
	replace p_em_coef = r(estimate) in `=`t'+2'
	replace p_em_se = r(se) in `=`t'+2'
	replace p_em_pvalue = 2 * ttail(r(df),abs(r(estimate)/r(se))) in `=`t'+2'
	
	local for_lincom `for_lincom' + _b[l`=`t'+1'.precip_CRU_pw50] + (2 * _b[l`=`t'+1'.precip_CRU_pw502] * 9)
}

* Lincom LIC
local for_lincom _b[precip_CRU_pw50] + (2 * _b[precip_CRU_pw502] * 11)
forval t = 0/$k {
	lincom `for_lincom'
	replace p_lic_coef = r(estimate) in `=`t'+2'
	replace p_lic_se = r(se) in `=`t'+2'
	replace p_lic_pvalue = 2 * ttail(r(df),abs(r(estimate)/r(se))) in `=`t'+2'
	
	local for_lincom `for_lincom' + _b[l`=`t'+1'.precip_CRU_pw50] + (2 * _b[l`=`t'+1'.precip_CRU_pw502] * 11)
}


export excel period - r2 in 1/9 using "output/Table_1.xlsx", sheet("column_8") sheetreplace firstrow(var)

